home *** CD-ROM | disk | FTP | other *** search
/ FM Towns: Free Software Collection 10 / FM Towns Free Software Collection 10.iso / ms_dos / tool / mls / mls.for < prev    next >
Text File  |  1995-02-17  |  20KB  |  642 lines
  1. ***********************************************************************
  2. ********************                        ***************************
  3. **************             MLS Ver.1.00          **********************
  4. ********************                        ***************************
  5. ***********************************************************************
  6. C
  7.       PROGRAM MLS
  8.       PARAMETER (MAX=100)
  9.       DIMENSION A(MAX,MAX),B(MAX,MAX),C(MAX,MAX),R(MAX,MAX),W(MAX,MAX)
  10.       DIMENSION D(MAX,MAX),EN(MAX),Z(MAX),XX(MAX),XS(MAX),V(MAX)
  11.       INTEGER X(MAX),JCH(MAX)
  12.       CHARACTER QT*1,FILE*25
  13.  1    ID=0
  14.       CALL START(MENU)
  15.       IF(MENU.EQ.5) THEN
  16.          CALL MAKE2(B,XX,K,N,D,JCH,MAX,MENU)
  17.          CALL COPYMTX(B,A,N,K,MAX)
  18.          CALL COPYMTX(XX,XS,N,1,MAX)
  19.          GOTO 46
  20.  111  END IF
  21.       IF(MENU.NE.2) THEN
  22.          CALL MAKE(A,W,XS,M,N,MENU,MAX,ICH)
  23.       ELSE
  24.  18      WRITE(6,104)'DATA 数を入力して下さい'
  25.          READ(5,*,ERR=18) M
  26.          IF(M.LE.0) THEN
  27.             WRITE(6,*) ' 入力しなおしてください！'
  28.             GOTO 18
  29.  24      END IF
  30.          WRITE(6,*)'DATA の入力方法を指定してください'
  31.  73      WRITE(6,*)'1:直接入力  2:DATA FILEによる入力'
  32.          READ(5,*,ERR=73) IDATA
  33.          IF(IDATA.LE.0.OR.IDATA.GE.3) GOTO 73
  34.  74      IF(IDATA.EQ.2) THEN
  35.             WRITE(6,*)'FILE の名前を入力してください'
  36.             READ(5,301) FILE
  37.             OPEN(1,FILE=FILE)
  38.          ELSE
  39.             WRITE(6,*)' DATA を X,Y の順に入力して下さい'
  40.          END IF
  41.          DO 4 I=1,M
  42.          IF(IDATA.EQ.1) THEN
  43.                READ(5,*) XX(I),EN(I)
  44.          ELSE
  45.                READ(1,*) XX(I),EN(I)
  46.                IF(I.EQ.M) CLOSE(1)
  47.          END IF
  48.  4       CONTINUE
  49.          WRITE(6,105)'反復数、傾きの順に出力します'
  50.          CALL BEW(XX,EN,XS,V,Z,M)
  51.          GOTO 68
  52.  5       CONTINUE
  53.       END IF
  54.       IDA=0
  55.       K=0
  56.       IF(ICH.EQ.1) THEN
  57.          K=M-1
  58.       END IF
  59.  11   K=K+1
  60.  12   CALL PRDMTX(W,A,B,N,N,K,MAX)
  61.       CALL PRDMTX(W,XS,XX,N,N,1,MAX)
  62.  46   CALL QR(B,R,X,EN,N,K,MAX,ID)
  63.       IF(ID.EQ.1) GOTO 7
  64.  25   CALL TRANSMTX(B,C,X,N,K,MAX)
  65.       CALL TENMTX(C,B,N,K,MAX)
  66.       CALL PRDMTX(B,XX,Z,K,N,1,MAX)
  67.       CALL BACK(R,XX,Z,K,MAX)
  68.       IF(MENU.EQ.1.AND.K.EQ.M) GOTO 6
  69.       IF(IDA.EQ.1.OR.ICH.EQ.1) GOTO 6
  70.       IF(MENU.EQ.5.OR.MENU.EQ.6) GOTO 6
  71.       CALL PRDMTX(W,A,C,N,N,K,MAX)
  72.       CALL TRANS(XX,V,X,K)
  73.       CALL COPYMTX(V,XX,K,1,MAX)
  74.       CALL PRDMTX(W,XS,V,N,N,1,MAX)
  75.       CALL COPYMTX(V,XS,N,1,MAX)
  76.       CALL PRDMTX(C,XX,V,N,K,1,MAX)
  77.  28   CALL AIC(V,XS,K,N,AI)
  78.       IF(MENU.NE.1) THEN
  79.          IF(K.EQ.1) THEN
  80.             IAI=K
  81.             SAI=AI
  82.          ELSE IF(AI+1.0.LT.SAI) THEN
  83.             IAI=K
  84.             SAI=AI
  85.          END IF
  86.       END IF
  87.       IF(K.NE.M) GOTO 11
  88.  17   K=IAI
  89.       IDA=1
  90.       IF(MENU.EQ.3) THEN
  91.          L=K
  92.       ELSE IF(MENU.EQ.4) THEN
  93.          L=K-1
  94.       END IF
  95.       IF(ICH.EQ.2) THEN
  96.          WRITE(6,201) L
  97.       END IF
  98.       GOTO 12
  99.  6    CALL SOLVE(XX,X,K,MENU)
  100.       IF(MENU.EQ.5) THEN
  101.  13      CALL TRANS(XX,V,X,K)
  102.          CALL COPYMTX(V,XX,K,1,MAX)
  103.          CALL PRDMTX(A,XX,V,N,K,1,MAX)
  104.          CALL AIC5(D,XS,XX,JCH,K,N,MAX,AI)
  105.       END IF
  106.  68   WRITE(6,106)'以上の結果になりました'
  107.  7    WRITE(6,103)'再び実行しますか？ ( Y or N )'
  108.       READ(5,108) QT
  109.       IF(QT.EQ.'N'.OR.QT.EQ.'n') THEN
  110.          CALL MES
  111.          STOP
  112.       ELSE IF(QT.EQ.'Y'.OR.QT.EQ.'y') THEN
  113.          GOTO 1
  114.  8    ELSE
  115.          GOTO 7
  116.  9    END IF
  117.       STOP
  118.  103  FORMAT(2X,A34)
  119.  104  FORMAT(/,2X,A27)
  120.  105  FORMAT(/,2X,A32)
  121.  106  FORMAT(/,5X,A35)
  122.  108  FORMAT(A1)
  123.  201  FORMAT( /,10X,I3,' 次式が最も良い近似です')
  124.  301  FORMAT(A25)
  125.       END
  126. ************************************************************************
  127. *                  SUBROUTINE TO OUTPUT THE MESSAGES                   *
  128. ************************************************************************
  129.       SUBROUTINE MES
  130.          WRITE(6,*)'THANK YOU FOR USING MLS!    SEE YOU AGAIN!!'
  131.          WRITE(6,*)'MERCI BEAUCOUP!!          AU REVOIR!!'
  132.          WRITE(6,*)'DANKE!!                   AUF WIEDERSEHEN!!'
  133.          WRITE(6,*)'謝謝!!                    再見!!'
  134.          WRITE(6,*)'GRAZIE!!                  ARRIVEDERCI!!'
  135.          WRITE(6,*)'GRACIAS!!                 HASTA LUEGO!!'
  136.          WRITE(6,*)'OBRIGADO!!                ATE LOGO!!'
  137.          WRITE(6,*)'СПАСИБО!!'
  138.          WRITE(6,*)'БАЯРЛАЛАА!!'
  139.       RETURN
  140.       END
  141. ************************************************************************
  142. *                 SUBROUTINE TO CHOOSE THE PROGRAM                     *
  143. ************************************************************************
  144.       SUBROUTINE START(MENU)
  145.       CHARACTER QT*1
  146.       WRITE(6,101) '最小二乗法プログラム'
  147.  19   WRITE(6,102)'近似する方法を選んで下さい'
  148.       WRITE(6,*)
  149.      *'1:多項式近似  2:比例式近似  3:連立方程式  4:重回帰モデル'
  150.       READ(5,*,ERR=19) MENU
  151.       IF(MENU.LE.0.OR.MENU.GE.5) GOTO 19
  152.       IF(MENU.EQ.4) THEN
  153.           MENU=5
  154.           RETURN
  155.       END IF
  156.       MENU=4-MENU
  157.       IF(MENU.EQ.3) THEN
  158.  2       WRITE(6,103)'定数項は入れますか？( Y or N )'
  159.          READ(5,108) QT
  160.          IF(QT.EQ.'Y'.OR.QT.EQ.'y') THEN
  161.             MENU=MENU+1
  162.          ELSE IF(QT.NE.'N'.AND.QT.NE.'n') THEN
  163.             GOTO 2
  164.  3       END IF
  165.       END IF
  166.       RETURN
  167.  101  FORMAT(/,5X,A46)
  168.  102  FORMAT(/,2X,A30)
  169.  103  FORMAT(2X,A34)
  170.  108  FORMAT(A1)
  171.       END
  172. ***********************************************************************
  173. *      SUBROUTINE TO IMPUT THE DATA AND TO MAKE THE MATRIX            *
  174. ***********************************************************************
  175.       SUBROUTINE MAKE(A,W,X,M,N,MENU,MAX,ICH)
  176.       DIMENSION A(MAX,MAX),X(1),W(MAX,MAX)
  177.       CHARACTER FILE*25
  178.       IF(MENU.EQ.1) THEN
  179.  11      WRITE(6,*)'パラメーターの数を入力して下さい'
  180.          READ(5,*,ERR=11) M
  181.          IF(M.LE.0) THEN
  182.             WRITE(6,*)' 入力しなおしてください！'
  183.             GOTO 11
  184.  34      END IF
  185.          N=M
  186.          WRITE(6,*)'各列毎に、左辺の係数、右辺の数の順に入力して下さい'
  187.          DO 2 I=1,M
  188.             DO 1 J=1,M
  189.                W(J,J)=1.0
  190.                READ(5,*) A(I,J)
  191.  1          CONTINUE
  192.             READ(5,*) X(I)
  193.  2       CONTINUE
  194.       ELSE
  195.  12      CONTINUE
  196.          WRITE(6,*) ' DATA 数を入力して下さい'
  197.          READ(5,*,ERR=12) N
  198.          IF(N.LE.0) THEN
  199.             WRITE(6,*) ' 入力しなおしてください！'
  200.             GOTO 12
  201.  29      END IF
  202.  30      WRITE(6,*)' 最高次数決定の方法を選択して下さい'
  203.          WRITE(6,*)' 1:直接入力  2:自動決定'
  204.          READ(5,*,ERR=30) ICH
  205.          IF(ICH.LE.0.OR.ICH.GE.3) THEN
  206.             WRITE(6,*) ' 選択し直して下さい'
  207.             GOTO 30
  208.  31      END IF
  209.          IF(ICH.EQ.2) THEN
  210.             M=INT(2*SQRT(REAL(N)))
  211.          ELSE
  212.  9          WRITE(6,*) ' 求める最高次数を入力して下さい'
  213.             READ(5,*,ERR=9) M
  214.             IF(M.GE.N) THEN
  215.                WRITE(6,*)'次数を入力しなおして下さい'
  216.                GOTO 9
  217.  7          END IF
  218.             IF(MENU.EQ.4) THEN
  219.                M=M+1
  220.             END IF
  221.          END IF
  222.          WRITE(6,*)'DATA の入力方法を指定してください'
  223.  23      WRITE(6,*)'1:直接入力  2:DATA FILEによる入力'
  224.          READ(5,*,ERR=23) IDATA
  225.          IF(IDATA.LE.0.OR.IDATA.GE.3) GOTO 23
  226.  24      IF(IDATA.EQ.2) THEN
  227.             WRITE(6,*)'FILE の名前を入力してください'
  228.             READ(5,101) FILE
  229.             OPEN(1,FILE=FILE)
  230.          ELSE
  231.             WRITE(6,*)' X,Y の順に入力して下さい'
  232.          END IF
  233.          DO 5 J=1,N
  234.             W(J,J)=1.0
  235.             IF(IDATA.EQ.1) THEN
  236.                READ(5,*) XDATA,YDATA
  237.             ELSE
  238.                READ(1,*) XDATA,YDATA
  239.             END IF
  240.             IF(MENU.EQ.4) THEN
  241.                DO 3 I=1,M
  242.                   A(J,I)=XDATA**(I-1)
  243.  3             CONTINUE
  244.             ELSE
  245.                DO 4 I=1,M
  246.                   A(J,I)=XDATA**I
  247.  4             CONTINUE
  248.             END IF
  249.             X(J)=YDATA
  250.  5       CONTINUE
  251.       END IF
  252.       IF(IDATA.EQ.2) CLOSE(1)
  253.       RETURN
  254.  101  FORMAT(A25)
  255.       END
  256. ***********************************************************************
  257. *                    SUBROUTINE TO EXECUTE QR-DIVISION                *
  258. ***********************************************************************
  259.       SUBROUTINE QR(A,R,X,EN,NA,AM,MAX,ID)
  260. C     MODIFIED GRAM-SCHMIDT METHOD
  261.       INTEGER AM
  262.       DIMENSION A(MAX,MAX),R(MAX,MAX),EN(1)
  263.       INTEGER X(1)
  264.       DO 1 I=1,NA
  265.       X(I)=I
  266.       DO 1 J=1,AM
  267.       R(I,J)=0.0
  268.  1    CONTINUE
  269.       DO 2 J=1,AM
  270.          CALL NORM(A,EN,X,J,NA,AM,MAXNORM,MAX)
  271.          IF(J.NE.AM) THEN
  272.             ISTORE=X(J)
  273.             X(J)=MAXNORM
  274.             IX=X(MAXNORM)
  275.             X(IX)=ISTORE
  276.             IF(J.GE.2) THEN
  277.                IF(X(J).NE.X(IX)) THEN
  278.                   DO 7 I=1,J-1
  279.                      STORE=R(I,X(X(J)))
  280.                      R(I,X(X(J)))=R(I,X(X(IX)))
  281.                      R(I,X(X(IX)))=STORE
  282.  7                CONTINUE
  283.                END IF
  284.             END IF
  285.          END IF
  286.          R(J,J)=EN(X(J))
  287.          IF(R(J,J).LT.1.0E-6) THEN
  288.              WRITE(*,*)'ランク落ちがおこりました'
  289.              WRITE(*,*)'もう一度やり直してください'
  290.              ID=1
  291.              RETURN
  292.          END IF
  293.          DO 3 I=1,NA
  294.             A(I,X(J))=A(I,X(J))/R(J,J)
  295.  3       CONTINUE
  296.          DO 4 K=J+1,AM
  297.             DO 5 L=1,NA
  298.                R(J,K)=R(J,K)+A(L,X(J))*A(L,X(K))
  299.  5          CONTINUE
  300.             DO 6 L=1,NA
  301.                A(L,X(K))=A(L,X(K))-R(J,K)*A(L,X(J))
  302.  6          CONTINUE
  303.  4       CONTINUE
  304.  2    CONTINUE
  305.       RETURN
  306.       END
  307. ***********************************************************************
  308. *            SUBROUTINE TO CALCULATE THE EUCLID NORM                  *
  309. ***********************************************************************
  310.       SUBROUTINE NORM(A,EN,X,IS,MA,NA,MAXNORM,MAX)
  311.       DIMENSION A(MAX,MAX),EN(1)
  312.       INTEGER X(1)
  313.       DO 1 JJ=IS,NA
  314.          J=X(JJ)
  315.          EN(J)=0.0
  316.          DO 2 I=1,MA
  317.             EN(J)=EN(J)+A(I,J)*A(I,J)
  318.  2       CONTINUE
  319.          EN(J)=SQRT(EN(J))
  320.          IF(JJ.EQ.IS) THEN
  321.             MAXNORM=J
  322.          ELSE
  323.             IF(EN(J).GT.EN(MAXNORM)) THEN
  324.                MAXNORM=J
  325.             END IF
  326.          END IF
  327.  1    CONTINUE
  328.       RETURN
  329.       END
  330. ***********************************************************************
  331. *                SUBROUTINE TO MAKE TRANSPOSED-MATRIX                 *
  332. ***********************************************************************
  333.       SUBROUTINE TENMTX(A,B,M,N,MAX)
  334.       DIMENSION A(MAX,MAX),B(MAX,MAX)
  335.       DO 1 I=1,N
  336.          DO 2 J=1,M
  337.             B(I,J)=A(J,I)
  338.  2       CONTINUE
  339.  1    CONTINUE
  340.       RETURN
  341.       END
  342. ***********************************************************************
  343. *        SUBROUTINE TO CALCCULATE THE PRODUCTION OF MATRIX A,B        *
  344. ***********************************************************************
  345.       SUBROUTINE PRDMTX(A,B,C,L,M,N,MAX)
  346.       DIMENSION A(MAX,MAX),B(MAX,MAX),C(MAX,MAX)
  347.       DO 30 I=1,L
  348.          DO 20 J=1,N
  349.             C(I,J)=0.0
  350.             DO 10 K=1,M
  351.                C(I,J)=C(I,J)+A(I,K)*B(K,J)
  352.  10         CONTINUE
  353.  20      CONTINUE
  354.  30   CONTINUE
  355.       RETURN
  356.       END
  357. ***********************************************************************
  358. *               SUBROUTINE TO TRANSLATE THE MATRIX                    *
  359. ***********************************************************************
  360.       SUBROUTINE TRANSMTX(A,B,X,M,N,MAX)
  361.       DIMENSION A(MAX,MAX),B(MAX,MAX)
  362.       INTEGER X(1)
  363.       DO 1 I=1,M
  364.          DO 2 J=1,N
  365.             B(I,J)=A(I,X(J))
  366.  2       CONTINUE
  367.  1    CONTINUE
  368.       RETURN
  369.       END
  370. ***********************************************************************
  371. *          SUBROUTINE TO EXECUTE THE BACKWARD SUBSTITUTION            *
  372. ***********************************************************************
  373.       SUBROUTINE BACK(R,A,Z,M,MAX)
  374.       DIMENSION R(MAX,MAX),A(1),Z(1)
  375.       DO 10 I=1,M
  376.          A(I)=0.0
  377.  10   CONTINUE
  378.       A(M)=Z(M)/R(M,M)
  379.       DO 1 I=M-1,1,-1
  380.          DO 2 J=I+1,M
  381.             A(I)=A(I)+R(I,J)*A(J)
  382.  2       CONTINUE
  383.          A(I)=(Z(I)-A(I))/R(I,I)
  384.  1    CONTINUE
  385.       RETURN
  386.       END
  387. ***********************************************************************
  388. *                    SUBROUTINE TO COPY MATRIX                        *
  389. ***********************************************************************
  390.       SUBROUTINE COPYMTX(A,B,M,N,MAX)
  391.       DIMENSION A(MAX,MAX),B(MAX,MAX)
  392.       DO 1 I=1,M
  393.          DO 2 J=1,N
  394.             B(I,J)=A(I,J)
  395.  2       CONTINUE
  396.  1    CONTINUE
  397.       RETURN
  398.       END
  399. ***********************************************************************
  400. *                      SUBROUTINE TO SOLVE THE MLS                    *
  401. ***********************************************************************
  402.       SUBROUTINE SOLVE(A,X,M,MENU)
  403.       DIMENSION A(1)
  404.       INTEGER X(1)
  405.       WRITE(*,203)
  406.       WRITE(*,*)
  407.       DO 10 LOOP=1,M
  408.          DO 20 J=1,M
  409.             IF(X(J).EQ.LOOP) THEN
  410.                I=J
  411.                GOTO 30
  412.  15         END IF
  413.  20      CONTINUE
  414.  30      CONTINUE
  415.          IF(MENU.EQ.1) THEN
  416.             WRITE(6,200) LOOP,A(I)
  417.          ELSE IF(MENU.EQ.3) THEN
  418.             WRITE(6,201) LOOP,A(I)
  419.          ELSE
  420.             WRITE(6,202) LOOP-1,A(I)
  421.          END IF
  422.  200     FORMAT('           X(',I3,') =',F13.3)
  423.  201     FORMAT('          ',I3,'次の係数 =',F13.6)
  424.  202     FORMAT('          ',I3,'次の係数 =',F13.6)
  425.  203     FORMAT( /,10X,'*****  結 果 の 出 力  *****')
  426.  10   CONTINUE
  427.       RETURN
  428.       END
  429. ***********************************************************************
  430. *                   SUBROUTINE TO CALCULATE MENU 2                   *
  431. ***********************************************************************
  432.       SUBROUTINE BEW(X,Y,W,SUM,Z,M)
  433.       DIMENSION X(1),Y(1),Z(1),W(1),SUM(1)
  434.       DO 1 I=1,M
  435.          W(I)=1.0
  436.          SUM(I)=W(I)
  437.  1    CONTINUE
  438.       DO 2 ICOUNT=1,10
  439.          XY=0.0
  440.          X2=0.0
  441.          DO 3 I=1,M
  442.             XY=XY+W(I)*X(I)*Y(I)
  443.             X2=X2+W(I)*X(I)*X(I)
  444.  3       CONTINUE
  445.          V2=0.0
  446.          A=XY/X2
  447.          WRITE(*,101) ICOUNT,A
  448.          DO 7 I=1,M
  449.             WRITE(6,102) I,W(I)
  450.  7       CONTINUE
  451.          DO 4 I=1,M
  452.             W(I)=ABS(Y(I)-A*X(I))
  453.             V2=V2+W(I)*W(I)
  454.  4       CONTINUE
  455.          CALL MEDIAN(W,M,Z,S)
  456.          C=20.8-15.0*EXP(-1.3730E-3*(S**1.5))
  457.          DO 5 I=1,M
  458.             IF(W(I).LT.S*C) THEN
  459.                W(I)=(1-(W(I)/(C*S))**2)**2
  460.             ELSE
  461.                W(I)=0.0
  462.             END IF
  463.  5       CONTINUE
  464.          HS=0.0
  465.          HB=0.0
  466.          DO 6 I=1,M
  467.             HS=HS+ABS(W(I)-SUM(I))
  468.             HB=HB+W(I)
  469.             SUM(I)=W(I)
  470.  6       CONTINUE
  471.          IF(HS.LT.HB*1.0E-4) RETURN
  472.  2    CONTINUE
  473.       RETURN
  474.  101  FORMAT(/,10X,I4,'回目の傾きの値 ＝',F13.6)
  475.  102  FORMAT(14X,I4,'番目の重み ＝',F13.6)
  476.       END
  477. ***********************************************************************
  478. *               SUBROUTINE TO SEARCH FOR THE MEDIAN                   *
  479. ***********************************************************************
  480.       SUBROUTINE MEDIAN(W,M,X,S)
  481.       DIMENSION W(1),X(1)
  482.       DO 1 I=1,M
  483.          X(I)=I
  484.  1    CONTINUE
  485.       ID=MOD(M,2)
  486.       IS=(M-ID)/2+1
  487.       DO 2 I=1,IS
  488.          DO 3 J=1,M-I
  489.             IF(W(X(J)).LT.W(X(J+1))) THEN
  490.                LARGE=X(J)
  491.                X(J)=X(J+1)
  492.                X(J+1)=LARGE
  493.             END IF
  494.  3       CONTINUE
  495.  2    CONTINUE
  496.       IF(ID.EQ.1) THEN
  497.          S=W(X(IS))
  498.       ELSE
  499.          S=(W(X(IS-1))+W(X(IS)))/2.0
  500.       END IF
  501.       RETURN
  502.       END
  503. ************************************************************************
  504. *              SUBROUTINE TO CALCULATE THE NUMBER OF AIC               *
  505. ************************************************************************
  506.       SUBROUTINE AIC(V,XS,M,N,AI)
  507.       DIMENSION XS(1),V(1)
  508.       PAI=3.145927
  509.       SV=0.0
  510.       DO 1 I=1,N
  511.          V(I)=XS(I)-V(I)
  512.          SV=SV+V(I)*V(I)
  513.  1    CONTINUE
  514.       SV=SV/REAL(N)
  515.       IF(SV.EQ.0.0) THEN
  516.          WRITE(*,*)
  517.          WRITE(*,*) '   残差０です！'
  518.          AI=-1.0E7
  519.          GOTO 5
  520.  2    END IF
  521.       AI=REAL(N)*(LOG(2*PAI*SV)+1.0)+2*(REAL(M)+1.0)
  522.   5   RETURN
  523.       END
  524. ************************************************************************
  525. *        SUBROUTINE TO TRANSLATE THE MATRIX ( SINGLE DIMENSION )       *
  526. ************************************************************************
  527.       SUBROUTINE TRANS(A,V,X,M)
  528.       DIMENSION A(1),V(1)
  529.       INTEGER X(1)
  530.       DO 1 LOOP=1,M
  531.          DO 20 J=1,M
  532.             IF(X(J).EQ.LOOP) THEN
  533.                I=J
  534.                GOTO 30
  535.  15         END IF
  536.  20      CONTINUE
  537.  30      V(LOOP)=A(I)
  538.  1    CONTINUE
  539.       RETURN
  540.       END
  541. ************************************************************************
  542. *                    SUBROUTINE TO EXECUTE MENU 5                      *
  543. ************************************************************************
  544.       SUBROUTINE MAKE2(X,Y,M,N,XDATA,JCH,MAX,MENU)
  545.       DIMENSION X(MAX,MAX),XDATA(MAX,MAX),Y(1)
  546.       INTEGER JCH(1)
  547.       CHARACTER FILE*25
  548.  1    WRITE(6,*) ' DATA 数を入力してください'
  549.       READ(5,*,ERR=1) N
  550.       IF(N.LE.0) GOTO 1
  551.          WRITE(6,*) ' パラメータ数を入力してください'
  552.          READ(5,*) M
  553.       IST=M
  554.  2    WRITE(6,*) ' DATA 入力の方法を入力してください'
  555.  3    WRITE(6,*) ' 1:直接入力  2:DATA FILE から入力'
  556.       READ(5,*,ERR=3) ICH
  557.       IF(ICH.LE.0.OR.ICH.GE.3) GOTO 3
  558.  4    IF(ICH.EQ.2) THEN
  559.         WRITE(6,*) ' FILE 名を入力してください'
  560.         READ(5,101) FILE
  561.         OPEN(1,FILE=FILE)
  562.       END IF
  563.       IF(MENU.EQ.5) THEN
  564.          JID=0
  565. C         IF(ICH.EQ.2) THEN
  566.          WRITE(6,*) '回帰するパラメータに１をしないものに０を'
  567.          WRITE(6,*) '入力して下さい'
  568.          DO 7 I=1,IST
  569.  6          WRITE(6,102) I
  570.             READ(5,*) JCH(I)
  571.             IF(JCH(I).NE.0.AND.JCH(I).NE.1) GOTO 6
  572.             IF(JCH(I).EQ.1) JID=JID+1
  573.  7       CONTINUE
  574.          M=JID+1
  575. C         ELSE
  576. C         DO 21 I=1,IST
  577. C             JCH(I)=1
  578. C 21      CONTINUE
  579. C         END IF
  580.          DO 8 I=1,N
  581.             IF(ICH.EQ.2) THEN
  582.                READ(1,*) (XDATA(I,J),J=1,IST),Y(I)
  583.             ELSE
  584.                WRITE(6,*) 'X,Y の順に入力してください'
  585.                READ(5,*) (XDATA(I,J),J=1,IST),Y(I)
  586.             END IF
  587.  8        CONTINUE
  588.           IF(ICH.EQ.2) CLOSE(1)
  589.           DO 22 I=1,N
  590.              X(I,1)=1.0
  591.  22       CONTINUE
  592.           ICOUNT=1
  593.           DO 11 J=1,N
  594.              DO 10 I=2,M
  595.  9              CONTINUE
  596.                 IF(JCH(ICOUNT).EQ.0) THEN
  597.                    ICOUNT=ICOUNT+1
  598.                    GOTO 9
  599.                 END IF
  600.                 X(J,I)=XDATA(J,ICOUNT)
  601.                 ICOUNT=ICOUNT+1
  602.  10          CONTINUE
  603.           ICOUNT=1
  604.  11       CONTINUE
  605.        END IF
  606.        RETURN
  607.  101   FORMAT(A25)
  608.  102   FORMAT(' X(',I3,')')
  609.        END
  610. ************************************************************************
  611. *          SUBROUTINE TO CALUCULATE THE NUMBER OF AIC AT MENU 5        *
  612. ************************************************************************
  613.       SUBROUTINE AIC5(X,Y,A,JCH,M,N,MAX,AI)
  614.       DIMENSION X(MAX,MAX),Y(1),A(1)
  615.       INTEGER JCH(1)
  616.       PAI=3.141593
  617.       SIGONE=0.0
  618.       SIGTWO=0.0
  619.       SIGTHR=0.0
  620.       DO 1 I=1,N
  621.          SIGONE=SIGONE+Y(I)**2
  622.          SIGTWO=SIGTWO+Y(I)
  623.  1    CONTINUE
  624.       SIGTWO=A(1)*SIGTWO
  625.       ICOUNT=1
  626.       DO 3 J=2,M
  627.          DO 2 I=1,N
  628.  10         IF(JCH(ICOUNT).EQ.0) THEN
  629.                ICOUNT=ICOUNT+1
  630.                GOTO 10
  631.  11         END IF
  632.             SIGTHR=SIGTHR+X(I,ICOUNT)*Y(I)*A(J)
  633.  2       CONTINUE
  634.          ICOUNT=ICOUNT+1
  635.  3    CONTINUE
  636.       SIGMA=(SIGONE-SIGTWO-SIGTHR)/REAL(N)
  637.       AI=N*(LOG(2*PAI)+1)+N*LOG(SIGMA)+2*(M+1)
  638.       WRITE(6,100) AI
  639.       RETURN
  640.  100  FORMAT( /,'            A I C =',F13.6)
  641.       END
  642.